Examples and Plotting
Plots, simpleยค
Plots.jl and Makie.jl are fully supported by Rasters.jl, with recipes for plotting Raster and RasterStack provided. plot will plot a heatmap with axes matching dimension values. If mappedcrs is used, converted values will be shown on axes instead of the underlying crs values. contourf will similarly plot a filled contour plot.
Pixel resolution is limited to allow loading very large files quickly. max_res specifies the maximum pixel resolution to show on the longest axis of the array. It can be set manually to change the resolution (e.g. for large or high-quality plots):
using Rasters, RasterDataSources, ArchGDAL, Plots
A = Raster(WorldClim{BioClim}, 5)
plot(A; max_res=3000)
For Makie, plot functions in a similar way. plot will only accept two-dimensional rasters. You can invoke contour, contourf, heatmap, surface or any Makie plotting function which supports surface-like data on a 2D raster.
To obtain tiled plots for 3D rasters and RasterStacks, use the function Rasters.rplot([gridposition], raster; kw_args...). This is an unexported function, since we're not sure how the API will change going forward.
Makie, simpleยค
using CairoMakie
CairoMakie.activate!(px_per_unit = 2)
using Rasters, CairoMakie, RasterDataSources, ArchGDAL
A = Raster(WorldClim{BioClim}, 5)
Makie.plot(A)

Loading dataยค
Our first example simply loads a file from disk and plots it.
This netcdf file only has one layer, if it has more we could use RasterStack instead.
using Rasters, NCDatasets, Plots
url = "https://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc";
filename = download(url, "tos_O1_2001-2002.nc");
A = Raster(filename)
180ร170ร24 Raster{Union{Missing, Float32},3} tos with dimensions:
X Mapped{Float64} Float64[1.0, 3.0, โฆ, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
Y Mapped{Float64} Float64[-79.5, -78.5, โฆ, 88.5, 89.5] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
Ti Sampled{DateTime360Day} DateTime360Day[DateTime360Day(2001-01-16T00:00:00), โฆ, DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals
extent: Extent(X = (-0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))missingval: missingcrs: EPSG:4326
mappedcrs: EPSG:4326
parent:
[:, :, 1]
-79.5 -78.5 โฆ 86.5 87.5 88.5 89.5
1.0 missing missing 271.43 271.437 271.445 271.459
3.0 missing missing 271.431 271.438 271.445 271.459
5.0 missing missing 271.431 271.438 271.445 271.459
7.0 missing missing 271.431 271.439 271.446 271.459
โฎ โฑ โฎ
351.0 missing missing 271.43 271.435 271.445 271.459
353.0 missing missing 271.43 271.436 271.445 271.459
355.0 missing missing 271.43 271.436 271.445 271.459
357.0 missing missing 271.43 271.437 271.445 271.459
359.0 missing missing โฆ 271.43 271.437 271.445 271.459
[and 23 more slices...]
Objects with Dimensions other than X and Y will produce multi-pane plots. Here we plot every third month in the first year in one plot:
A[Ti=1:3:12] |> plot
Now plot the ocean temperatures around the Americas in the first month of 2001. Notice we are using lat/lon coordinates and date/time instead of regular indices. The time dimension uses DateTime360Day, so we need to load CFTime.jl to index it with Near.
using CFTime
A[Ti(Near(DateTime360Day(2001, 01, 17))), Y(-60.0 .. 90.0), X(45.0 .. 190.0)] |> plot
Now get the mean over the timespan, then save it to disk, and plot it as a filled contour.
Other plot functions and sliced objects that have only one X/Y/Z dimension fall back to generic DimensionalData.jl plotting, which will still correctly label plot axes.
using Statistics
# Take the mean
mean_tos = mean(A; dims=Ti)
180ร170ร1 Raster{Union{Missing, Float32},3} tos with dimensions:
X Mapped{Float64} Float64[1.0, 3.0, โฆ, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
Y Mapped{Float64} Float64[-79.5, -78.5, โฆ, 88.5, 89.5] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
Ti Sampled{DateTime360Day} DateTime360Day(2002-01-16T00:00:00):Millisecond(2592000000):DateTime360Day(2002-01-16T00:00:00) ForwardOrdered Explicit Intervals
extent: Extent(X = (-0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))missingval: missingcrs: EPSG:4326
mappedcrs: EPSG:4326
parent:
[:, :, 1]
-79.5 -78.5 โฆ 86.5 87.5 88.5 89.5
1.0 missing missing 271.427 271.434 271.443 271.454
3.0 missing missing 271.427 271.434 271.443 271.454
5.0 missing missing 271.427 271.434 271.443 271.454
7.0 missing missing 271.427 271.435 271.444 271.454
โฎ โฑ โฎ
351.0 missing missing 271.427 271.433 271.443 271.454
353.0 missing missing 271.427 271.433 271.443 271.454
355.0 missing missing 271.427 271.433 271.443 271.454
357.0 missing missing 271.427 271.433 271.443 271.454
359.0 missing missing โฆ 271.427 271.433 271.443 271.454
Plot a contour plotยค
using Plots
Plots.contourf(mean_tos; dpi=300, size=(800, 400))
write to diskยค
Write the mean values to disk
write("mean_tos.nc", mean_tos)
"mean_tos.nc"
Plotting recipes in DimensionalData.jl are the fallback for Rasters.jl when the object doesn't have 2 X/Y/Z dimensions, or a non-spatial plot command is used. So (as a random example) we could plot a transect of ocean surface temperature at 20 degree latitude :
A[Y(Near(20.0)), Ti(1)] |> plot
Polygon masking, mosaic and plotยค
In this example we will mask the Scandinavian countries with border polygons, then mosaic together to make a single plot.
First, get the country boundary shape files using GADM.jl.
using Rasters, RasterDataSources, ArchGDAL, Shapefile, Plots, Dates, Downloads, NCDatasets
WARNING: using Downloads.download in module Main conflicts with an existing identifier.
Download the shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "boundary_lines.shp"
Downloads.download(shapefile_url, shapefile_name)
"boundary_lines.shp"
Load using Shapefile.jl
shapes = Shapefile.Handle(shapefile_name)
denmark_border = shapes.shapes[71]
norway_border = shapes.shapes[53]
sweden_border = shapes.shapes[54]
Polygon(4665 Points)
Then load raster data. We load some worldclim layers using RasterDataSources via Rasters.jl:
climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)
RasterStack with dimensions:
X Projected{Float64} LinRange{Float64}(-180.0, 179.833, 2160) ForwardOrdered Regular Intervals crs: WellKnownText,
Y Projected{Float64} LinRange{Float64}(89.8333, -90.0, 1080) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
:tmin Float32 dims: X, Y (2160ร1080)
:tmax Float32 dims: X, Y (2160ร1080)
:prec Int16 dims: X, Y (2160ร1080)
:wind Float32 dims: X, Y (2160ร1080)
mask Denmark, Norway and Sweden from the global dataset using their border polygon, then trim the missing values. We pad trim with a 10 pixel margin.
mask_trim(climate, poly) = trim(mask(climate; with=poly); pad=10)
denmark = mask_trim(climate, denmark_border)
norway = mask_trim(climate, norway_border)
sweden = mask_trim(climate, sweden_border)
RasterStack with dimensions:
X Projected{Float64} LinRange{Float64}(9.5, 25.6667, 98) ForwardOrdered Regular Intervals crs: WellKnownText,
Y Projected{Float64} LinRange{Float64}(70.5, 53.6667, 102) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
:tmin Float32 dims: X, Y (98ร102)
:tmax Float32 dims: X, Y (98ร102)
:prec Int16 dims: X, Y (98ร102)
:wind Float32 dims: X, Y (98ร102)
Then load raster data. We load some worldclim layers using RasterDataSources via Rasters.jl:
climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)
RasterStack with dimensions:
X Projected{Float64} LinRange{Float64}(-180.0, 179.833, 2160) ForwardOrdered Regular Intervals crs: WellKnownText,
Y Projected{Float64} LinRange{Float64}(89.8333, -90.0, 1080) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
:tmin Float32 dims: X, Y (2160ร1080)
:tmax Float32 dims: X, Y (2160ร1080)
:prec Int16 dims: X, Y (2160ร1080)
:wind Float32 dims: X, Y (2160ร1080)
mask Denmark, Norway and Sweden from the global dataset using their border polygon, then trim the missing values. We pad trim with a 10 pixel margin.
mask_trim(climate, poly) = trim(mask(climate; with=poly); pad=10)
denmark = mask_trim(climate, denmark_border)
norway = mask_trim(climate, norway_border)
sweden = mask_trim(climate, sweden_border)
RasterStack with dimensions:
X Projected{Float64} LinRange{Float64}(9.5, 25.6667, 98) ForwardOrdered Regular Intervals crs: WellKnownText,
Y Projected{Float64} LinRange{Float64}(70.5, 53.6667, 102) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
:tmin Float32 dims: X, Y (98ร102)
:tmax Float32 dims: X, Y (98ร102)
:prec Int16 dims: X, Y (98ร102)
:wind Float32 dims: X, Y (98ร102)
Plotting with Plotsยค
First define a function to add borders to all subplots.
function borders!(p, poly)
for i in 1:length(p)
Plots.plot!(p, poly; subplot=i, fillalpha=0, linewidth=0.6)
end
return p
end
borders! (generic function with 1 method)
Now we can plot the individual countries.
dp = plot(denmark)
borders!(dp, denmark_border)
and sweden
sp = plot(sweden)
borders!(sp, sweden_border)
and norway
np = plot(norway)
borders!(np, norway_border)
The Norway shape includes a lot of islands. Lets crop them out using .. intervals:
norway_region = climate[X(0..40), Y(55..73)]
plot(norway_region)
And mask it with the border again:
norway = mask_trim(norway_region, norway_border)
np = plot(norway)
borders!(np, norway_border)
Now we can combine the countries into a single raster using mosaic. first will take the first value if/when there is an overlap.
scandinavia = mosaic(first, denmark, norway, sweden)
RasterStack with dimensions:
X Projected{Float64} 3.1666666666666443:0.16666666666666666:32.49999999999998 ForwardOrdered Regular Intervals crs: WellKnownText,
Y Projected{Float64} 72.66666666666666:-0.16666666666666666:52.99999999999999 ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
:tmin Float32 dims: X, Y (177ร119)
:tmax Float32 dims: X, Y (177ร119)
:prec Int16 dims: X, Y (177ร119)
:wind Float32 dims: X, Y (177ร119)
And plot scandinavia, with all borders included:
p = plot(scandinavia)
borders!(p, denmark_border)
borders!(p, norway_border)
borders!(p, sweden_border)
p
And save to netcdf - a single multi-layered file, and tif, which will write a file for each stack layer.
write("scandinavia.nc", scandinavia)
write("scandinavia.tif", scandinavia)
(tmin = "scandinavia_tmin.tif", tmax = "scandinavia_tmax.tif", prec = "scandinavia_prec.tif", wind = "scandinavia_wind.tif")
Rasters.jl provides a range of other methods that are being added to over time. Where applicable these methods read and write lazily to and from disk-based arrays of common raster file types. These methods also work for entire RasterStacks and RasterSeries using the same syntax.
This page was generated using Literate.jl.